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ABSTRACT 


Computer Simulations were performed to locate the equilibrium 
positions and binding energies of interstitial He, Ne, Ar, Kr, and 
Xe atoms ina tungsten crystal. Heavy interstitial atoms in tung- 
sten share a lattice site with the atom that normally occupies that 
E. and form what is called a split interstitial. Three charac- 
teristic interstitial sites were located relative to each lattice 
site tested. The distance of the impurity atom from the site was 
Seen to vary roughly inversely with its mass, and the dues CONDE 
of the lattice atom increased with the mass of the impurity atom. 

The foreign atom in its interstitial position was tested to 
determine the minimum initial kinetic energy needed to escape the 
lattice, as well as the optimum escape direction. The minimum 
energy may be interpreted to be the binding energy of the defect. 
A comparison of experimental binding energies from Kornelsen and 
Sinha and simulated binding energies indicates the model gives 


realistic results. 
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I. INTRODUCTION 


The advent of modern, high speed digital computers has led to 
the application of computer simulation techniques to many differ- 
ent types of physical systems. One such application is the modeling 
of the situation that occurs when a foreign atom interacts with a 
metallic crystal lattice. In general, such modeling can be broken 
down into two basic areas, dynamic simulation and static simulation. 
AS an example of the former, radiation damage has been studied by 
the simulated firing of an atom or ion onto a crystal face. Other 
examples are sputtering simulations [1,2,3] in which the incoming 
particle causes surface atoms to be ejected; and channeling simu- 
lations [4] in which the ranges of ions travelling in crystal 
lattices are calculated. Static simulations, on the other hand, 
have been concerned with the equilibrium BEE in the lattice 
after point defects such as replacement atoms, interstitial atoms, 
and vacancies have been introduced. Examples of this type of 
simulation can be found in απ. present research uti- 
lized aspects of both static and dynamic simulation techniques. 

The goal of this research was to correlate the results of experi- 
mentally determined binding energies of point defects in a tungsten 


lattice [8,9] with the results obtained by computer simulation. 





II. NATURE OF THE PROBLEM 


A. HISTORICAL BACKGROUND 

An investigation of radiation damage events by computer simu- 
lation techniques of crystalline behavior was published in 1960 by 
Gibson, Goland, Milgram, and Vineyard, hereafter referred to as 
(GGMV) [5]. This Brookhaven National Laboratory investigation set 
forth the criteria that must be satisfied before the simulation 
method could be applied to crystals. Such factors as potential 
energy functions, forces, crystallite sizes, computation methods, 
choice of time intervals, and computer limitations were discussed. 
The crystal lattice modeled in their research was metallic copper, 
a face-center cubic (fcc) structure. The potential functions used 
in the calculations was a Born-Mayer repulsive potential, with the 
necessary cohesive forces applied on the boundries of the crystal- 
lite. In integrating the equations of motion, the Brookhaven 
group used the central difference procedure. The optimum choice 
for timestep duration, At, was shown to be of critical importance 
in the integration scheme. The energy of the strongest interac- 
tion governed their choice of the above parameter. The static 
results obtained by GGMV confirmed the existance of the (100) 
split interstitial in the fcc lattice. Their dynamic results 
described collision chains and focusing phenomena in crystallites 
Struck with energetic knock-on atoms. 

Additional crystal simulation studies were performed by 


R.A. Johnson ο το 11). In Ref. 6, he investigated point defects 


in a copper lattice using Born-Mayer repulsive potentials. The 
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(100) split interstitial was found to be the only stable inter- 
stitial position. He found it necessary to allow the interstitial 
to interact only with its six nearest neighbors. Atoms near the 
defect were treated as independent, while the remainder of the 
metal was treated as an elastic continuum with atoms imbedded in 
те. 

Research on body-centered cubic (bcc) crystals was undertaken 
by Erginsoy, Vineyard, and Englert (EVE) IL... They used a com- 
posite potential function for most of their detailed work. It 
consisted of an exponentially screened Coulomb potential at small 
separations, a Born-Mayer function in the region between small 
and intermediate separations, and a Morse potential at larger 
separation distances. A split interstitial was reported for simu- 
lated bee crystals in the (110) direction. In their dynamic re- 
sults they reported the existence of a threshold energy for 
displacement that was highly independent of the direction of 
knock-on. Also, collisional chains in the (111) and (100) direc- 
tions were found to exist. R.A. Johnson also published results 
for bcc simulations involving @-iron and tungsten Ree 
existence of the split interstitial in the (110) direction was 
confirmed and crowdian migration data were discussed. The present 
investigation confirmed the (110) split interstitial positions 
for argon, neon, krypton, and xenon. 

D.E. Harrison and associates have published several articles 
in which computer simulation of crystalline behavior has been 
investigated Шоо]; In a study of a fcc model of copper, col- 


lision events between a copper atom and a copper lattice were 


10 
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simulated as a function of the potential function, the enegy of 
the collision, and the location of the impact point. The inte- 
gration scheme used was an average force procedure, instead of 
the central difference procedure used by GGMV. A complete dis- 
cussion of the method can be found in Ref. 13. A continuation of 
copper simulations was published by Harrison, Leeds, and Gay in 
1965 i». Another paper by Harrison and Greiling dealt with the 
ranges of heavy ions in tungsten crystals whose atoms had under- 
gone thermal displacement (ს). It was found that room temperature 
thermal displacements had a negligible effect upon the collisions 
for ion energies greater than a few thousand electron volts. 
Finally, Harrison, Levy, Johnson, and Effron published results on 
computer simulation of sputtering [2]. 

Research undertaken by Vine [14] provided the foundation for 
this author's present investigation. Vine used the Gay-Harrison 
model for crystal simulation as modified by Levy [15], Johnson [16], 
Effron [17], and Moore [18]. His overall objective was to cor- 
relate experimental and simulated binding energies of neon and 
argon M ibt defects. Repulsive potential functions for neon- 
tungsten (Ne-W) and argon-tungsten (Ar -W) interactions were 
used [19]. Morse functions were not used for those interactions 
because experimental data giving the Morse parameters was based on 
homogeneous media such as tungsten-tungsten (W-W) interactions [20]. 
Vine subsequently attempted to correlate results of equilibrium 
position studies for tungsten defects in a tungsten lattice using 
two different potential functions for the interstitial-lattice 


interactions. In one case, the tungsten interstitial was allowed 
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to interact using a Born-Mayer repulsive potential, and the other 
case, the tungsten defect was given a composite potential that was 
identical to that given to all the other lattice atoms. By so 
doing, he attempted to establish an empirical relationship between 
the two methods to apply to the results of neon interstitial studies 
which used only a repulsive potential. 

The results of the tungsten-tungsten interactions failed to 
provide the information needed to formulate a correction factor 
to be used in the neon-tungsten studies. In addition, the concept 
of relating the potential energy at equilibrium to the experi- 
mentally observed binding energies of Kornelsen and Sinha [9] does 


not appear to be feasible. 


B. CHOICE OF POTENTIAL FUNCTION 

The studies that were reviewed in the previous section utilized 
many different approximations to the true potential function be- 
tween atoms in a metal lattice. The problem of solving the many- 
body interaction of a real system is approximated in the computer 
by many two-body interactions. Thus, central pairwise potential 
functions are most often used in computer simulations. GGMV 


employed a repulsive potential of the Born-Mayer form: 


"as = oe eh 


which described the repulsion of atoms at close approach. Three 
Born-Mayer potentials were investigated by GGMV to determine which 
would give the best results in their calculations. Their choice 

was one which has since been labeled the Gibson Number Two Potential. 


In Ref. 6, Johnson and Brown used a similar Born-Mayer potential 
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with slightly different parameters. As previously mentioned, 
another potential function that has been used in crystal modeling 


studies is the familiar Morse potential of the form: 


Var - D| exp (- 20%, , - D -2 ep (af, - 291 


where rg is the equilibrium distance of approach of two atoms, 


and Œ and D are constants. 

Girifalco and Weizer calculated Morse parameters that would 
be appropriate for several crystal lattices o calculating 
the parameters, they attempted to express the various crystal 
properties such as cohesive energy, lattice constant, compres- 
sibility and equation of state in terms of the Morse function. 
The Morse potential constants published by Girifalco and Weizer 
have been used extensively in simulating the potential functions 
and forces in lattices of homonuclear atom systems. 

The Born-Mayer potential and the Morse potential are useful 
over specific internuclear separation distances. The Born-Mayer 
potential is useful at strongly repulsive separations, i.e. short 
ranges, while the Morse function is applicable at equilibrium and 
greater separations. In order to better approximate the true 
potential function, various investigators have combined different 
potential functions. As mentioned earlier, EVE [7] used a com- 
posite potential that consisted of an exponentially screened 
Coulomb potential at very close separations, a Born-Mayer po- 
tential at weakly repulsive distances, and a Morse potential over 
the remainder of the potential curve. In their studies, Harrison 


and associates combined a Born-Mayer potential and a Morse 
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potential. These two potentials were joined by a cubic equation 


in the region near their intersection В А 


С. EXPERIMENTAL RESULTS OF KORNELSEN AND SINHA 

In 1968, Kornelsen and Sinha [7,8] published results concerning 
the binding energies of trapped particles such as neon, argon, 
krypton, and xenon ions in a tungsten surface. The particles were 
forced into the surface from a beam created by an ion gun which 
gave ion energies of 40 eV to 5 keV. The tungsten crystal was 
then heated and the rates of evolution of the trapped gas were 
measured. 

The temperature at which desorption peaks occured thus gave 
an indication of the binding energies of point defects in the 
tungsten crystal.  Quantitatively similiar results were obtained 
with argon, krypton, neon, and xenon ions. Four peaks were ob- 
served below 1650 Әк апа were labeled as Q peaks. A single peak 
above 1700 °K was measured and was called the m It was con- 
cluded that the O-group of peaks were the results of single step 
desorptions from sites very close to the metal surface. They 
further concluded that the different &@-peaks could correspond to 
different types of point defect binding energies in the tungsten 
crystal. Specifically mentioned were defects of three types; 
(a) interstitial and substitutional positions in the lattice, (b) 
different distances from the site to the surface, and (c) different 
locations of a nearby lattice defect, such as a vacancy, relative 


to the site and the surface. 


14 





III. THE SIMULATION MODEL 


A. THE CRYSTAL 

The model used in this research was essentially the same one 
used by Vine (14) as explained in Section II. In subsequent dis- 
cussion, when it is necessary to specify certain of the computer 
program variable names, they will be placed in braces. 

All of the simulations in this research were done with tungsten 
crystals of varying sizes. The objective was to use the smallest 
crystal dimensions that would give realistic results. The di- 
mensions described below refer to the number of planes of atoms 
in each of the three rectangular coordinate directions. Simu- 
lations were performed with sizes 8 x 6 x 8 of 96 atoms, 

10 x 6 x 10 of 150 atoms, 10 x 8 x 10 of 200 atoms, and 

10 x 10 x 10 of 250 atoms. Of these, the latter two were judged 

to be of most use because of their greater depth in the y-direction. 
The y direction was always used as the direction of escape for the 
point mercer atom. 

Each atom in the simulated crystal was numbered, with the 
first position always assigned to the point defect atom. The 
remainder of the positions were assigned in sequence according to 
the coordinate locations. The numbering was started in the y = O 
plane and continued until all the atoms in that plane were speci- 
fied. This procedure was repeated for the remainder of the y planes 


in the simulated crystal. 


15 





B. THE TIMESTEP INTERVAL 

The numerical method of time integration used in the model was 
the average force method eae The value of At {pr} used in this 
procedure was of critical importance in determining whether or not 
the model would approximate reality. Also, the program running 
time was a function of the timestep duration. 

In order to best approximate reality, the {DT} was kept smaller 
early in the program when most movement was expected, and was al- 
lowed to grow larger as equilibrium was approached. The parameter 
that controlled the timestep duration was (ori): the EF cT E, 
the most energetic atom was allowed to move in a single timestep. 
In previous work, (DTI) was held constant throughout the duration 
of a program. If the motion was expected to be slow, (ртт) was 
given a larger value than in a situation where simulated motion 
was expected to be greater. For static equilibrium problems, at 
least 100 timesteps were needed to reach an approximately stable 
position. In addition, it became necessary to insert a maximum 
value for the timestep, поле into the program to prevent unrea- 
listic movement of the atoms and breakdown of the model. 

In the current program, changes were made to allow {ртт} to 
decrease during the program. For example, in static runs on He-W, 
Ne-W, 'Ar-W, and Kr -W, (отт) was initially given a value of 0.1 
lattice unit ee Tungsten forms a bcc crystal, with LU equal 
ВОС or 1.598 А, where LC is the Lattice Constant or cube edge 
distance. {DTI} was allowed to change in decrements of 0.01 LU 
per timestep for the first 10 timesteps. Then the {DTI} value 


of 0.01 LU was allowed to change in decrements of .001 LU for 
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another 10 timesteps. Equilibrium was reached by 30 timesteps 


according to this procedure, resulting in a significant decrease 


of computer running time. Also, the equilibrium positions obtained 


in the simulation were closer to the expected (110) split inter- 


stitial positions than were obtained with a constant {ртт} value. 
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ТУУ IMULATION PROCEDURE 


A. STATIC SIMULATIONS 

The first stage of the simulation procedure was concerned with 
finding the equilibrium positions for the point defect of interest 
in the top several layers of a bcc tungsten crystal. Previous work 
done by Johnson fill, indicated that a bcc crystal had at least 
two different types of stable split interstitial orientations. 
The first of these was in a (100) direction along a (110) plane. 
The other stable orientation was in a (110) direction in a (100) 
plane. An analysis of the total of 12 stable split interstitial 
positions possible about a given atom in the two orientations 
listed above, revealed that there were only three independent types 
of sites. (See Figure 1.) The first of these is located ona 
(110) plane with {Nvac} and lies closer to the surface than {nvac}. 
The second independent position lies in the same (100) plane as 
{ МУАС} and 1s at the same depth in the crystal. A final position 
similar to the first, in that it is along a (110), is located at 
a depth below that of (NVAC]. The three different split inter- 
stitials were postulated to have different binding energies be- 
cause of their varying depths in the crystal. 

In order to reduce the running time in locating the final 
coordinate positions at each interstitial location, the static 
programs were started with the foreign atoms in the approximate 


area of the expected final positions as, discussed above. For 
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example, when it was desired to locate position A, the foreign 
atom was initially placed along the correct (110) plane approxi- 


mately one lattice unit from {NvAc}. 


B. DYNAMIC SIMULATIONS 

The dynamic simulation program used the final positions com- 
puted in the static programs as input initial crystal positions. 
The lattice generator subroutine {B100}, and the point defect 
locator {PLACE} were thus eliminated from the dynamic program. 

Dynamic simulations were broken down into two main categories. 
The first category consisted of survey runs of the areas above 
the interstitials in the direction of escape from the y surface. 
An impact point generator package was inserted into the dynamic 
program, thus permitting the interstitial to be directed at a 
Specific number of locations in a predefined area. The results 
of the survey run were analyzed to determine the optimum aiming 
point for the interstitial at a specific initial energy. It was 
observed during the impact testing precedure that the "best" 
aiming points were a function of the initial energy given to the 
interstitial. However, it was also found that the optimum aiming 
points at varying energies were in a generally localized area. 
Thus, once the optimum points were found at the starting energy, 
only a localized region around those points was tested at lower 
initial interstitial energies. The decrementing process of the 
initial interstitial energy constituted the second main phase of 
the dynamic simulations. The procedure in this phase was to de- 


crease the interstitial energy until it could no longer escape 
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from the crystal. The minimum escape energy was said to be the 
simulated binding energy of the interstitial at the particular 


location tested. 
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V. PRESENTATION OF DATA 


OSTATIC RUNS 
Preliminary Testing of Static Program 

One hundred twelve computer runs were made in the initial 
testing phase of the static program. The PUIK of the testing was 
done in four general areas: 1) Program shutdown procedures at 
equilibrium; 2) Optimum crystal size determination; 3) Equilibrium 
position as a function of the initial interstitial position; and 
4) Realistic timestep determination procedure. 

The best method of stopping the program at equilibrium has 
been a matter of იარ for many years with the static simulation 
program. In the present investigation, the first method tested was 
a shutdown procedure initiated whenever a sharp {DT} decrease was 
encountered in the program. This test proved to be of limited 
success, and was later abandoned. Another method that was at- 
tempted was a test of { EMAX} against a value such as .04, to deter- 
mine if Ва eg had been reached. This method also proved only 
partially successful. Finally, a test was made of the average 
kinetic energy of an atom in the simulated crystal. The average 
energy was taken to be the total kinetic energy {TPKE} divided by 
the number of atoms ure or equivalently, ( TPKE) multiplied by 
the reciprocal of the number of atoms (RLL). The crystal was 
assumed to be at equilibrium if the average kinetic energy of an 
atom was less than or equal to the value 0.025 eV, a value for the 


average thermal energy of an atom. Satisfactory results were 
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sometimes obtained by this method but the test had to be removed 
in many cases to allow the simulated crystal to run a greater 
number of timesteps and reach equilibrium. 

As previously mentioned, different sizes of crystals were 
simulated to determine the optimum dimensions for the crystallite. 
Many tests were made on crystal sizes smaller than the 10 x 10 x 10 
used by Vine μα... Although smaller sizes such as 10 x 6 x 10 often 
gave reliable results, it was finally decided to use the 10 x 10 x 10 
size for the reported split interstitial positions in order to have 
a standard size applicable over all crystal positions of interest. 

In his work, Vine placed his interstitials in the middle 
of the open channels in the simulated crystal. Numerous runs in 
the present research have indicated that equilibrium is reached 
more quickly when the initial starting positions are not in 
channel centers, but along the directions of the expected splits 
in the general area of the final positions. 

All of the initial work was done with a fixed {DTI} value 
in the program. The { DTI} value normally used was .05 LU. The 
timestep Ee was later modified as explained in The Simulation 
Model and the computer running time was cut by approximately two- 
thirds due to this procedure. 

2. Positions for Helium in Tungsten 

Helium was the lightest point defect used in the static 
Simulation model. It was thought that the small, light atom would 
essentially do all the movement and come to rest in a channel 
center V 2 Lხ away from ( NVAC] along the (110? direction. The 


results of the runs show that the movement was almost as expected. 
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Also, a comparison of the C site for atom 64 and the A site for 
atom 114 indicates that there is a strong possibility that these 
two sites are degenerate. This is based on the fact that the 
potential energies of both sites are close, and also that the 
( NVAC] for each site is displaced only a negligible distance. The 
same probability of degeneracy is also seen to exist for C-89 and 
A-139. The determination of degeneracy of sites is seen to bea 
complex evaluation of potential energy differences, movements of 
{nvac}, and the potential gradient between the two sites. The 
final positions for the A, B, and C split interstitial positions 
in the third, fourth, fifth, and sixth layers where {NVAC} was 
64, 89, 114, and 139 respectively are shown in Figures 2, 3 and 4. 
The results obtained with helium were thought to be somewhat in 
error, since the heavier atom, neon, moved further in its simu- 
lations than helium. Further work is needed to confirm the equi- 
librium positions for helium. 
3. Positions for Neon in Tunasten 

The neon split interstitial locations were simulated for 
the same positions as helium. (See Figures 5, 6 and 7.) Essenti- 
ally, the ( NVAC) atom remained in its lattice site and the neon 
moved along a (110) direction to a distance N 2 LU away from 
(NvAC). 

ს. Positions for Argon in Tungsten 

The bulk of the testing of this research was done with 

argon as the simulated point defect. The results of the runs for 


the 10 x 10 x 10 crystal size are shown'in Figures 8, 9 and 10. 
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Figure 3 
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Figure 5 
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Figure 6 
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Figure 7 
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Figure 8 
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Figure 9 
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Figure 10 
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The expected (110) splits were observed with {nvac} moving away 


from its site in relation to the final.position of the argon defect. 


- 


5. Positions for Krypton in Tungsten 


The split interstitial positions for krypton ina simu- 
lated tungsten crystal were also plotted. (See Figures 11, 12 
and 13.) Once e Satisfactory results were obtained and the 
(110) splitting was observed. 

6. Positions for Xenon in Tungsten 

For the Xenon Runs, the tungsten repulsive potential 
function was used for the xenon potential. The initial runs, 
using the (pri) employed for the other defects, failed to give 
the postulated split interstitial positions. At least two factors 
were seen to complicate the Xenon-Tungsten runs. Firstly, the 
mass of Xenon was approximately 6/7 that of tungsten, so the 
lattice site would be shared almost evenly. This would mean that 
tungsten would have to move a significant distance from its lat- 
tice position. Secondly, the relatively large size of the xenon 
defect would require more movement of the surrounding atoms to 
accommodate the extra atom. 

The { ртт) scheme was changed to allow for more movement of 
the atoms by allowing (pri) to change in smaller decrements. 
Another method that was used was an initial displacement of both 
the xenon and the tungsten away from the lattice site. When both 
atoms were displaced, it became necessary to use an initial {ptr} 
of .05 LU or less. An average value of the C site for atom 89 
was computed. (See Figure 14.) More work is needed to simulate 


all of the interstitial positions. 
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Figure 12 


да C SITES FOR ATOMS 89, 139 


TUNGSTEN 29 | 





(Y)LU 40 





Figure 13 


| B SITES FOR ATOMS 64, 89, 114, 139 


KRYPTON-QO | TUNGSTEN — C) 





37 





Figure 14 
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opt interstitial Distance Ratios as a Function of 

Relative Masses 

During the static testing phase, the ratios of the split 
interstitial distances from the initial ( NVAC) Were measured. An 
attempt was made to correlate these ratios to the inverse ratios 
of the atomic masses of the atoms involved. (See Table I.) The 
point defect results of helium and neon did not give any signi- 
ficant correlation, but the Argon and krypton atoms gave ratios 
of splits in good agreement with the expected values from mass 


ratio calculations. 


B. DYNAMIC RUNS 


1. Survey of the Possible Directions of Escape for the Ar-W 
Simulated Split Interstitial in Site A-89 


As was mentioned in the Section IV, survey runs were made 
of the possible escape directions in the plane one LU. above the 
defect. The results for the survey for the Ar-W split interstitial 
in site A-89 were plotted. (See Figure 15.) {py}, the distance in 
lattice units that the argon moved in 25 timesteps, and {vy}, the 
velocity that the argon had after 25 timesteps are shown for each 
impact point. Due to the multiple scattering of the defect as it 
moved toward the surface, 25 timesteps were not sufficient to 
provide conclusive evidence at any one impact point. The survey 
was limited to 25 timesteps per impact point, however, due to 
computer time considerations, The main benefit gained from the 
survey run was the elimination of certain areas from further testing, 
Such as those on the outer perimeter of the survey area. Also, much 


information was gained on the most likely area of escape for the 
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TABLE I 
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Tungsten Distance 
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Figure 15 
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defect. From the results of the survey run, and from symmetry 
considerations of the open channel, it was decided to concentrate 
further testing on several points with the same Z coordinate as 
the interstitial aton. 
2. Determination of the Simulated Binding Energy of Argon 

in Site A-89 

Detailed testing of Several impact points with the same Z 
coordinate as the argon was carried out for 100 timesteps per 
impact point. This was generally enough time for’ the simulated 
>" лат to escape the crystal it its initial energy was 


great enough. The results are summarized below. 


TABLE II Detailed Impact Point Testing 





Impact Points in 


Initial Ener y =>” 






Plane 1.0 LU in 20 eV 10 ev S ev 6 ev ს CV 
-y direction from 
Yes Yes No No No 
Yes Yes Yes Yes No 
Yes Yes Yes Yes No 


Yes - Atom Escapes Crystal 


Νο - Atom Does Not Escape 


42 





VI. CONCLUSIONS 


The method used to decrease {ртт} on every timestep appears 
to be a Significant improvement over the old method of keeping 
{ртт} fixed. In effect, the new method forces the simulated 
crystallite to come to equilibrium in a predefined number of 
timesteps. In addition to a saving of computer time, this method 
also eliminates the need for a equilibrium shut down procedure 
in the program. The same {γττ) decrement scheme could not be used 
for all the point defects tested in this research. For example, 
Argon and Krypton were able to come to their expected equilibrium 
positions using the scheme described in Section III, but Xenon 
and Helium were not. Thus it appears that defect size, as well 
as the expected degree of movement may dictate the {DTI} decre- 
ment scheme. 

The expected (110) split interstitials were SE for the 
helium, neon, argon, krypton, and xenon foreign defects. The re- 
lative degree of splitting does appear to be related to the masses 
of the interacting atoms, but conclusive evidence of this was only 
obtained for the EECH and Krypton-Tungsten runs. 

The dynamic runs have shown the direction of escape from the 
crystal to be a multi-collisional process, with the defect under- 
going many intermediate direction changes at the low kinetic 
energies tested. The most likely escape direction was also seen 
to vary somewhat with the initial kinetic energy given to the 


defect. 
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The simulated value of approximately 4 eV for the binding 
energy of Argon in site A-89 appears to be a reasonable value. 
Kornelsen's data showed an energy level at this value and his 
levels were postulated to arise from defects in the first several 
layers of the tungsten crystal. 

Further work 1s needed to confirm the remainder of the energy 
levels found by Kornelsen, by testing other point defect sites in 


the top layers of the simulated crystal. 
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APPENDIX A: COMPUTER PROGRAM GLOSSARY 


NOTE: In this glossary, the terms "point defect atom", "bullet", 
and "primary" are synonymous; and the terms "latttice atom" and 


"target" are synonymous. 


AC: Distance measurement used in impact point generator 
ALPHA: Input Morse potential parameter 
BSAVE: Target mass/(target mass + bullet mass); distributed 


potential energy between target and bullet 


BIND: Negative of the total potential energy (TPOT) at time 
zero 
BMAS: Mass of bullet in amu 


BULLET: Alpha-numeric array for point defect material 

СЕО, СЕ|, CF2: Force parameters of cubic fit between Morse and 
Born-Mayer functions 

CGB1, CGB2: Morse potential parameters 

CGD1, CGD2: Morse potential parameters 

ССЕ1 , CGF2: Morse force parameters 

COX, COY, COZ: Cosines of angles to x,y, and z axes respectively 

СРО, СР1, CP2, CP3: Potential parameters of cubic fit between 
Morse and Вогп-Мауег functions 

CVD: CVR x 102, converts lattice units to meters 

CVE: 1.6 x ილა converts electron volts to joules 

CVED: CVE/CVD, a ratio used to avoid repeated division 

CVM: 1.672 x ara: converts atomic mass units to kilograms 


CVR: LU in angstroms; converts lattice units to angstrom units 
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ο. ο, ο Coordinates of impact point 
D1X, D1Y, D1Z: Displacement coordinates for location of inter- 


stitial from reference atom, NVAC 


DCON : Input Morse potential parameter 

DDII: Time increment that is subtracted from DTI after each 
timestep 

DFF: ROE-DIST, the distance closer than ROE that an atom is to 


the primary 

DIST: Distance between any two atoms 

DLPE: TLPE-TLPEØ, the change in total local potential energy 
since time zero 

DRX, DRY, DRZ: x,y,z components of DIST 

DT: Length of a timestep in seconds 

DTI: Number of lattice units most energetic atom may move in 


one timestep 


DITS: Starting value of DTI. 
DTOD: DT/CVD--a ratio used to avoid repeated division 
DTOM: DT/PTMAS--a ratio used to avoid repeated division 


DTOMB: DT/PEMAS--a ratio used to avoid repeated division 
DX(I), DY(I), DZ(I): Change in position of ith atom from initial 


position at time zero 


EMAX: The maximum energy encountered in any cycle 
ERAT: Measure of the average kinetic energy of an atom 
EV: Primary energy in electron volts 

EVR: Primary energy in kilo-electron volts 


EXA, EXB: Input Born-Mayer potential function parameters for the 


target 
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F2: 


FA: 


Square of the force on a specific atom 


The component force increment on an atom 


FDTI: DTI X CVD, a parameter used to determine DT by maximum 


FM2: 


FMAX: 


FOD: 


FORCE: 


FX(I), F 


FXA: 


HBMAS: 


HDTOD: 


HDTOM: 


HDTOMB: 


HTMAS: 


mI: 


ТЮБЕР: 


IH1: 


IH2: 


IHB: 


IHS: 


energy method 

A small number used in checking potential energy zero 
point 

FM squared 

Maximum total force on the most stressed atom in the 
crystal 

Maximum total force on Atom 1 

FORCE/DIST--a ratio used to avoid repeated division 
Numerical value of the force function with a variable 


parameter 


Y(I), FZ(1I): x,y,z components of total foce on an atom 


Born-Mayer force function parameter 


4s BMAS-- a ratio used to avoid repeated division 


4%; DIOD-- a ratio used to avoid repeated division 


% DTOM-- a ratio used to avoid repeated division 


- 


N 


% DTOMB-- a ratio used to avoid repeated division 


> TMAS-- a ratio used to avoid repeated division 
Variable in cubic fit subroutine 

Variable in cubic fit subroutine 

Number of mobile layers 

Alpha numeric array for program title 

Alpha numeric array for Morse function parameters 
Alpha numeric array for bullet element 


Alpha numeric array for type and orientation of crystal 
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ΠΤ: 


ELAY: 


IN: 


ВР = 


IO: 


ISHUT: 


EI: 


ETT: 


ETYPE: 


Alpha numeric array for target element 

Same as IDEEP 

Odd-even integer used to determine atom site establishment 
Subscript value of atom. Used in subroutine STEP and 

ENERGY 

Parameter that determines whether or not a self defect is 

to be given a repulsive potential or a composite attractive- 
repulsive potential 

A parameter used to shut down the program 

Unscaled fixed point x coordinate used in lattice generation 
Odd-even integer used to determine atom site establishment 
Parameter used to determine the type of point defect: 
vacancy, self-interstitial, replacement, foreign inter- 


Sti tial 


IVACX, IVACY, IVACZ: Input plane numbers to specify NVAC 


IX, IY, 
J2: 

JT: 
15: 
JTT: 
KF: 

KT: 


ЕСТ (І): 


LD: 


EI: 


IZ: Number of x,y,z planes of crystal 


Variable in the cubic fit subroutine 

Unscaled y coordinate used in crystal generation 
Variable used to establish atom sites 

Variable used to establish atom sites 

Final K in LOCAT (K) assigned to an atom 

Unscaled z coordinate used to establish atom site 
Used to identify an ith atom which is not included in 
calculations 


The highest numbered atom in the mobile layers 


The highest numbered atom in the entire crystal 
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LOCAT(K):Dimensioned variable that remembers the numbers of the 


MCRO: 


ND: 


NDEC: 


NEW: 


NPAGE: 


NRUN: 


PAC: 


PEMAS: 


atoms within a radius ROEL of the primary at time zero 
Variable associated with each of the nine lattice gene- 
rator subroutines 

One number higher than the order of the fit between the 
Born-Mayer and Morse potentials, always 4 in this simu- 
lation 

Data output increment, in numbers of timesteps 

Counting index for DTI variation 

Parameter used to determine whether or not atom numbers 
have been stored in LOCAT(K) 

Page numbering variable 

Parameter used to determine whether or not to read 
additional data cards 

Initial print statement timestep number 

Timestep number 

Timestep number limit before shutdown 

An atom number used to establish point defects or used as 
AN point for interstitial placement 

Parameter for bullet force function correction 


Primary mass in kilograms 


PEXA, PEXB: Input Born-Mayer potential function parameters for the 


PFPTC: 


PFXA: 


РКЕ(Т): 


PLANE: 


bullet-target interaction 

Primary force function evaluated at ROE 
Primary force function parameter 
Kinetic energy of the ith atom 


Alpha-numeric array for lattice orientation 
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POT: 


PPE(I): 


ЕРТС: 


PTE(I): 


PTMAS: 


RE: 


RLL: 


RO: 


ROE: 


ROE2: 


ROEA: 


ROEB: 


ROEC: 


КОЕС2 : 


ROEL: 


ROEL2: 


ROEM: 


Potential energy between two atoms 

Potential energy of the ith atom 

Primary potential function evaluated at ROE 

Total energy of the ith atom (potential + kinetic) 
Target mass in kilograms 

Input Morse potential parameter 

Reciprocal of LL 

Spacing constant in FCC(110) lattice generation subroutine 
Nearest neighbor distance 

ROE squared 

Maximum cut off for Born-Mayer potential 

Minimum cut off for Morse potential 

Maximum cut off for Morse potential 

ROEC squared 

Radius inside of which local potential energy is found 
ROEL squared 

ROE-DTI, region in which modification of repulsive force 


must be made 


RX(I), RY(I), RZ(I): x,y,z coordinates of an ith atom at any time 


ο ο τε ο RZI(I): x,y,Z coordinates of an ith atom's initial 


position 


RXK(I), RYK(I), RZK(I): x,y,z coordinates of temporary position of 


SAVE: 


an ith atom during force cycle 


12 POT 


SCX, SCY, SCZ: x,y,z coordinate scale factors 


START: 


SUM: 


An optional timing variable, not used in this simulation 


Variable in cubic fit subroutine 
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TARGET: 


TSAVE: 


TE: 


TEMP: 


TEFAC: 


TFAC: 


TFACB: 
THERM: 
TIME: 
TLPE: 
Т.РЕЯ: 
TMAS: 
TPKE: 
TPOT: 
TPOTL: 


VSS: 


Alpha-numeric array for target material 

Bullet mass/(target mass + bullet mass); distributes 
potential energy between target and bullet 

Total energy of all crystal atoms (kinetic + potential) 
Temperature of lattice in degrees Kelvin. Not used in 
this simulation 

Product of DTI * CVD 

A time factor ratio used to determine DI by maximum force 
method 

TFAC for the bullet 

Thermal energy of atom. Not used in this simulation 
Elapsed problem time in seconds 

Total local potential energy of atoms within a radius ROEL 
TLPE at time zero 

Target atom mass in amu 

Total kinetic energy of all crystal atoms 

Total potential energy of all crystal atoms 

Storage position for the last computed value TPOT 


Storage variable for velocity components 


VX(I), VY(I), VZ(I): x,y,z components of ith atoms velocity 


X, Y, Z: Unscaled coordinates used in crystal generation 


XSTART : 


X Coordinate used in impact point generator 


YLAX(I): Relaxation in -y direction of ith layer in L.U. 


ZE: 


Floating point form of JTT 
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